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Abstract 

We  present  two  stochastic  failure  models  for  the  reliability  evaluation  of  manufacturing 
equipment  that  degrades  due  to  its  complex  operating  environment.  The  first  model  examines 
the  case  when  the  environment  is  a  temporally  nonhomogeneous  continuous-time  Markov  chain, 
and  the  second  assumes  the  environment  is  a  temporally  homogeneous  semi-Markov  process  on 
a  finite  space.  Derived  are  transform  expressions  for  the  lifetime  distributions.  A  few  examples 
are  provided  to  illustrate  the  main  results. 
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1  Introduction 


Modern  manufacturing  organizations  employ  very  complex,  flexible  machines  that  can  be  used 
to  manufacture  a  variety  of  components  using  a  wide  range  of  material  types.  A  critical  factor 
in  the  economics  of  modern  manufacturing  facilities  is  the  useful  lifetime  of  machine  tools  that 
degrade  over  time  under  normal  usage  and  fluctuating  operating  conditions  that  can  be  driven  by 
external  influences.  It  is  well  known  that  many  factors  can  influence  the  rate  at  which  machine  tools 
degrade.  Some  of  the  most  important  factors  include  the  workpiece  material,  the  tool  material, 
the  cutting  speed,  and  the  depth  of  cut,  to  name  only  a  few.  Besides  these  obvious  contributors 
to  tool  degradation,  the  ambient  environment  of  the  manufacturing  equipment  can  significantly 
impact  the  rate  at  which  tools  (or  other  critical  components)  degrade.  For  instance,  significant 
fluctuations  in  the  ambient  temperature  or  humidity  can  cause  workpiece  (or  tool)  expansion  and 
contraction,  and  even  premature  corrosion.  Because  machine  tools  constitute  a  significant  portion 
of  a  manufacturing  organization’s  capital,  and  because  their  degradation  or  failure  can  significantly 
impact  production  and  inventory  levels,  it  is  important  to  assess  the  useful  lifetime  of  manufacturing 
equipment  when  it  is  subjected  to  complex  dynamics. 

In  this  paper  we  propose  two  stochastic  models  that  can  be  used  to  derive  mathematical  ex¬ 
pressions  for  the  reliability  of  components  that  degrade  dynamically  due  to  their  complex  ambient 
environment  or  time-varying  operating  conditions.  The  rate  of  degradation  is  modulated  by  an 
external  process  (the  random  environment)  which  is  modeled  as  a  continuous-time,  discrete-state 
stochastic  process  {Z(t)  :  t  >  0}.  At  some  initial  time,  say  s,  the  machine  tool  (or  component) 
begins  its  lifetime  in  perfect  working  order  but  degrades  over  time  under  the  influence  of  its  random 
environment  until  the  cumulative  degradation  reaches  a  fixed  threshold  x,  at  which  time  it  is  no 
longer  useful  because  it  may  not  be  able  to  maintain  manufacturing  specifications.  Derived  herein 
are  transform  expressions  for  the  cumulative  distribution  function  (c.d.f.)  of  the  useful  lifetime  of 
the  tool  when  the  evolution  of  the  operating  environment  assumes  two  distinct  forms.  First,  we 
consider  the  case  when  the  environment  evolves  as  a  nonhomogeneous,  continuous-time  Markov 
chain  (NHCTMC)  on  a  finite  state  space.  Second,  we  consider  the  case  of  a  time-homogeneous 
semi-Markov  process  (SMP)  environment,  also  on  a  finite  state  space.  For  the  first  model,  the  main 
transient  results  are  obtained  by  showing  that  the  joint  distribution  of  the  cumulative  degrada¬ 
tion  level  and  the  environment  state  conditionally  satisfies  a  nonhomogeneous  partial  differential 
equation.  This  equation  leads  to  the  Laplace-Stieltjes  transform  (LST)  of  the  lifetime  distribution 
function  under  certain  mild  conditions.  For  the  second  model,  we  employ  a  classical  Markov  renewal 
argument  to  derive  the  (double)  Laplace  transform  of  the  component’s  lifetime  distribution. 

Many  researchers  have  proposed  stochastic  models  of  machine  tool  reliability,  primarily  for 
the  purpose  of  prescribing  optimal  tool  replacement  policies,  or  to  minimize  processing  times.  A 
sampling  of  these  contributions  can  be  found  in  [10,  11,  12,  26,  30].  Others  have  considered  the 
control  of  either  manufacturing  or  production/inventory  systems  using  nonhomogeneous  Markov 
chains  to  model  the  failure  rate  which  is  assumed  to  be  dependent  on  the  production  rate  of  the 
system  (cf.  [24,  27,  28,  29]).  More  generally,  stochastic  failure  models  that  attempt  to  capture  the 
impact  of  a  randomly  evolving  environment  have  been  examined  extensively  in  the  reliability  and 
applied  probability  literature.  An  excellent  survey  due  to  Singpurwalla  [39]  presents  a  number  of 


2 


models  with  various  attributes.  Much  earlier  work,  due  to  Esary,  et  al.  [6]  provided  several  results 
for  wear  and  shock  processes.  Many  authors  have  considered  the  case  when  the  environment  evolves 
as  a  homogeneous,  continuous-time  Markov  chain.  For  example,  Qinlar  [4]  demonstrated  that  the 
joint  process  of  the  unit’s  wear  level  and  the  state  of  its  ambient  environment  may  be  considered 
as  a  Markov-additive  process.  He  considered  the  case  when  the  random  environment  is  a  general 
Markov  process  with  a  finite  state  space  and  the  wear  is  assumed  to  increase  as  a  Levy  process.  Li 
and  Luo  [23]  considered  a  Markov-modulated  shock  process  wherein  the  shock  inter-arrival  times 
and  the  random  shock  damage  are  both  governed  by  a  Markov  chain.  They  obtain  reliability  bounds 
when  the  inter-arrival  times  have  heavy-  or  light-tailed  distributions;  however,  their  degradation 
model  does  not  include  a  continuous  wear  component.  Kiessler  et  al.  [19]  investigated  the  limiting 
average  availability  of  a  system  whose  time-varying  wear  rates  are  governed  by  a  continuous-time 
Markov  chain.  Kharoufeh  et  al.  [15]  extended  the  model  of  [19]  by  including  damage-inducing 
shocks  that  arrive  according  to  a  time-homogeneous  Poisson  process  and  deriving  the  Laplace- 
Stieltjes  transforms  of  a  few  transient  reliability  indices.  Ebrahimi  [5]  considered  a  system  whose 
degradation  is  comprised  of  a  continuous  wear  component  as  well  as  jumps.  The  properties  of 
the  model  were  investigated  and  bounds  were  established  for  the  reliability  function.  However, 
the  random  environment  in  that  model  corresponds  to  a  shock  process  that  is  superimposed  on  a 
gamma  wear  process. 

Models  that  consider  temporally  nonhomogeneous  Markov  or  homogeneous  semi-Markov  envi¬ 
ronments  can  also  be  found  in  the  stochastic  modeling  and  operations  research  literature.  Platis  et 
al.  [35,  37,  38]  have  made  extensive  use  of  nonhomogeneous  Markov  models  for  assessing  the  perfor¬ 
mance  and  reliability  (or  performability)  of  complex,  time-dependent  systems,  while  Smotherman 
and  Zemoudeh  [40]  examined  phased-mission  reliability  models  using  nonhomogeneous  Markov 
models  to  reflect  globally-time-dependent  phase  changes.  Platis  et  al.  [36]  analyzed  hitting  times 
for  a  discrete-time  nonhomogeneous  Markov  chain  and  showed  how  these  can  be  used  for  reliability 
evaluation. 

The  analysis  we  present  here  serves  to  extend  the  results  of  [13,  14,  19]  and  others  by  considering 
environments  that  are  either  (i)  temporally  nonhomogeneous  in  their  stochastic  evolution,  or  (ii) 
non-Markovian  but  temporally  homogeneous.  Here,  we  solely  focus  our  attention  on  transient 
reliability  indices  (i.e. ,  when  the  critical  degradation  threshold  and  time  are  both  finite).  For 
the  NHCTMC  model,  we  provide  a  full  characterization  of  the  Laplace-Stieltjes  transform  of  the 
lifetime  distribution  when  the  infinitesimal  generator  matrix  entries  satisfy  some  mild  conditions. 
For  the  SMP  model,  we  provide  a  double  transform  characterization  of  the  lifetime  distribution 
and  illustrate  the  results  assuming  a  few  state  sojourn  time  distributions. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2  provides  a  brief  model  description 
followed  by  the  full  analysis  of  the  NHCTMC  model  in  section  3.  Section  4  presents  the  results  for 
the  SMP  model,  while  section  5  provides  a  few  concluding  remarks  and  some  directions  for  future 
research. 
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2  Model  Description 


Here,  we  describe  a  general  framework  for  modeling  equipment  degradation  in  the  context 
of  complex  manufacturing  environments.  In  the  following  discussion,  we  assume  that  all  random 
variables  are  defined  on  a  common  and  complete  probability  space  denoted  by  (P,  A,  P).  We  assume 
that  the  equipment  (hereafter  referred  to  as  the  unit)  is  placed  into  service  at  some  initial  time 
s  (s  >  0)  with  no  degradation  (i.e. ,  its  condition  is  new  initially).  Due  to  normal  usage  and 
the  influence  of  its  environment,  it  accumulates  degradation  until  a  deterministic,  critical  threshold 
value  x  is  reached  or  exceeded,  at  which  time  the  system  is  considered  to  be  failed.  As  a  prototypical 
example,  consider  manufacturing  operations  involving  metal  cutting  and  material  removal.  Most 
cutting  tools  degrade  in  a  manner  that  is  dictated  by  the  cutting  speed,  the  workpiece  material, 
cutting  angle,  use  of  cutting  fluids,  and  a  number  of  other  critical  factors.  However,  once  the  wear 
level  of  the  tool  reaches  a  critical  level  x,  it  may  not  be  possible  for  it  to  maintain  manufacturing 
specifications.  Our  primary  concern  is  the  development  of  a  stochastic  model  for  characterizing 
the  random  amount  of  time  needed  for  the  unit’s  cumulative  degradation  level  to  first  reach  this 
critical  threshold.  We  denote  this  time  by  a  proper  non- negative  random  variable  Tx(s).  We  include 
the  dependence  of  this  first  passage  time  on  the  initial  time  s;  however,  when  the  environment  is 
homogeneous  and  we  assume  s  =  0,  we  will  suppress  the  variable  s. 

The  degradation  of  the  unit  is  induced  by  (i)  the  actual  ambient  environment  in  which  it  resides 
and  operates  (e.g.,  the  ambient  air  temperature  or  humidity  can  induce  degradation);  and/or  (ii) 
the  various  operational  settings  of  the  equipment.  Degradation  of  the  first  type  might  occur  due 
to  temperature  fluctuations  that  cause  contraction  and  expansion  of  the  workpiece  and/or  the 
cutting  tool.  For  the  second  type,  if  different  materials  are  processed  according  to  an  a  priori 
schedule  (e.g.,  according  to  a  cyclic  Markov  chain),  but  the  processing  times  are  stochastic,  the 
unit  will  experience  different  rates  of  degradation  depending  on  the  current  workpiece.  The  key 
point  here  is  that  the  rate  of  degradation,  and  the  time  spent  in  each  degradation  rate,  are  governed 
by  a  continuous-time  stochastic  process  termed  the  random  environment  which  shall  be  denoted 
throughout  by  Z  =  {Z(t)  :  t  >  0}.  The  state  space  of  this  process  is  a  finite  set  S  =  {1,  2, . . .  ,£} 
where  each  state  of  S  corresponds  to  the  current  prevailing  environment  or  operating  condition. 

The  primary  objective  of  this  paper  is  to  analyze  two  specific  and  distinct  models  for  the 
governing  environment  Z.  For  the  first  model,  we  assume  that  Z  is  a  nonhomogeneous  continuous¬ 
time  Markov  chain  (NHCTMC)  whose  transition  functions  (and  infinitesimal  generator  matrix) 
depend  explicitly  on  time.  In  the  second  case,  we  will  examine  a  relaxation  of  the  Markovian 
assumption  by  allowing  Z  to  evolve  as  a  semi-Markov  environment  but  insist  that  the  environment 
be  time-homogeneous  for  this  latter  case.  For  both  models,  whenever  Z{t)  =  j  €  S,  the  unit 
degrades  at  a  constant  rate  r(j)  >  0.  More  specifically,  the  rate  of  degradation  evolves  as  a 
stochastic  process  and  is  modulated  by  the  environment.  For  later  use,  we  define  the  diagonal  matrix 
containing  these  environment-dependent  rates  by  =  diag(r(l),  r(2), . . .  ,r(£)).  The  assumption 
of  a  constant  degradation  rate  in  each  environment  state  amounts  to  assuming  that  the  degradation 
process  has  samples  paths  that  are  piecewise  linear  and  monotonically  increasing  almost  surely 
(a.s.).  This  assumption  is  not  restrictive  since  the  degradation  sample  path  can  be  approximated 
by  a  piecewise  linear  path  as  noted  in  [14].  Furthermore,  the  approach  is  well-suited  for  modeling 
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linear  degradation  paths  which  are  common  in  many  applications. 

For  both  the  NHCTMC  and  SMP  models,  given  that  the  unit  is  placed  into  service  at  time  s 
in  perfect  condition,  we  define  the  cumulative  degradation  of  the  unit  up  to  time  t  by  a  random 
variable  X(t).  For  any  t  >  s,  this  random  variable  is  given  by 

X(t)  =  X(s)  +  J  r(Z(v))dv  (1) 

where  we  assume  X(s)  =  0  a.s.  and 

r(Z(v))\dv  <  oo  a.s. 

so  that  X(t)  is  well  defined  for  each  t  >  0.  Because  we  assume  r(j)  >  0  for  each  j  €  S,  the  sample 
paths  of  { X(t )  :  t  >  0}  are  monotonically  increasing  a.s.,  and  consequently,  the  events  {X(t)  <  a:} 
and  {Tx(s)  >  t}  are  equivalent.  The  system’s  random  lifetime  is  the  first  passage  time 

Tx(s)  =  inf{i  >  0  :  X(t)  >  x},  (2) 

or  the  first  time  the  degradation  process  {X(i)  :  t  >  0}  reaches  or  exceeds  x.  Let  F(x,s,t )  = 
P (Tx(s)  <  t)  =  1  —  ¥(X(t)  <  x)  denote  the  unconditional  c.d.f.  of  the  unit’s  lifetime,  and  let  its 
conditional  counterpart  be  Fi(x,s,t)  =  P (Tx(s)  <  t\Z(s)  =  i),  i  G  S ,  where  i  denotes  the  state 
of  the  environment  at  the  initial  time  s.  In  the  case  when  the  environment  is  time-homogeneous, 
we  drop  the  dependence  on  s  and  simply  write  Tx,  F(x,t),  and  Fi(x,t )  for  the  first  passage  time, 
its  unconditional  c.d.f.,  and  its  conditional  c.d.f.,  respectively.  The  indicator  function  of  an  event 
A  G  A  will  be  denoted  by  1(A),  and  it  assumes  the  value  1  if  A  holds  and  0  otherwise.  Our 
primary  aim  is  to  derive  transform  expressions  for  the  unconditional  and/or  conditional  lifetime 
distributions  for  the  NHCTMC  and  SMP  models. 

3  Time-Nonhomogeneous  Markov  Environment 

In  this  section,  we  consider  the  case  in  which  the  unit’s  environment  evolves  as  a  finite,  non- 
homogeneous  continuous-time  Markov  chain  (NHCTMC).  These  processes  afford  a  great  deal  of 
modeling  flexibility  and  have  been  applied  in  a  variety  of  contexts  including  manufacturing  and 
production  systems  (cf.  [27,  28])  and  remaining  lifetime  prognosis  in  health  care  ([34]).  Here,  we 
focus  our  attention  on  transient  (as  opposed  to  asymptotic)  reliability  indices  in  order  to  account 
for  the  environment’s  explicit  dependence  on  time.  Our  first  main  result  shows  that  the  case  of 
a  time- homogeneous,  Markovian  environment  (cf.  [14,  15])  can  be  extended  very  naturally  to  the 
time-nonhonrogeneous  case. 

3.1  Transient  Lifetime  Distribution 

For  each  t  >  0,  let  Z(t)  denote  the  status  of  the  environment  at  time  t,  and  assume  Z  =  {Z(t)  : 
t  >  0}  is  a  NHCTMC  with  state  space  £  =  { 1,2,  E  Z+.  For  any  s,t  such  that  0  <  s  <  t, 

the  transition  functions  of  Z  are  given  by 

Pij(s,t)  =  F(Z(t)  =  j\Z{s )  =  i),  i,j  <E  S, 
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and  the  time-nonhomogeneous  transition  matrix  is  denoted  by  P(s,t )  =  \pij(s,  Note  that 

P(s,s )  =  I,  where  I  denotes  the  identity  matrix.  It  is  well  known  that  the  time-dependent  in¬ 
finitesimal  generator  matrix  of  Z,  denoted  by  Q(s),  s  >  0,  has  elements 

Qij(s)  =  lim  +  k\  i,j  G  S,  j  i, 

and  when  j  =  i, 

/  \  /  n  r  1  -Pii(s,s  +  h) 

=  Qi{s)  =  lim - - - 

h— >0  h 

so  that  qi(s)  =  These  entries  are  interpreted  as  follows:  When  j  i.  qij(s)  is  the 

instantaneous  rate  at  which  the  environment  transitions  from  state  i  to  state  j  at  time  s.  When 
j  =  i,  qi(s )  is  the  total  rate  at  which  the  process  transitions  out  of  state  i  at  time  s.  (Note 
that  qa{s)  <  0  for  all  s  >  0.)  In  case  the  environment  is  time-homogeneous,  then  of  course, 
P(s,  t )  =  P(t  —  s )  and  Q(s)  =  Q  for  all  s  >  0. 

For  any  s  and  t  such  that  0  <  s  <  t,  define 


Vij(x,s,t)  =  P(X(f)  <  x,Z(t)  =  j\Z(s)  =  i)  i,j  G  S, 


the  joint  probability  that,  at  time  t,  the  degradation  of  the  system  has  not  exceeded  level  x,  and 
the  environment  process  is  in  state  j  G  S ,  given  that  at  that  initial  time  s  the  environment  was  in 
state  i  G  S.  For  ease  of  notation,  we  also  define  the  t  x  t  matrix  V ( x ,  s,  t )  =  [Vij(x.  s,  t)].  As  noted 
in  section  2,  because  the  degradation  rates,  { r(j )  :  j  G  S},  are  assumed  to  be  strictly  positive  real 
numbers,  {X(t)  :  t  >  s}  is  a  continuous  additive  functional  of  Z  whose  sample  paths  are  monotone 
increasing  almost  surely.  This  monotonicity  ensures  that  {X(t)  <  x}  =  { Tx(s )  >  t};  therefore,  we 
can  derive  the  c.d.f.  (or  complementary  c.d.f.)  of  Tx(s)  by  analyzing  the  cumulative  degradation 
X(t).  To  this  end,  we  next  characterize  the  joint  distribution  functions  Vij(x,  s,t).  This  analysis 
is  similar  to  the  arguments  presented  in  [13,  14,  15],  but  requires  some  extra  care  to  account  for 
the  time-nonhomogeneity  of  Z.  The  following  theorem  shows  that  the  matrix  V(x,s,t )  verifies  a 
Chapman-Kolmogorov-type  (partial)  differential  equation. 


Theorem  1  If  the  operating  environment,  {Z(t)  :  t  >  0},  is  a  nonhomogeneous  CTMC  with  in¬ 
finitesimal  generator  matrix  Q(t),  then  the  matrix  V(x,  s,t )  satisfies  the  partial  differential  equation 


dV ( x ,  s,  t ) 

Ft 


,  dV(x,s,t)^ 

H - ^ - 1'/ 

ox 


V (x,  s,  t)Q{t). 


where  =  diag(r(l), . . . ,  r(£)). 


(3) 


Proof.  The  proof  of  the  result 
of  the  time-homogeneous  case  presented 


is  standard  for  Markov  processes  and  mirrors  the  proof 
in  [14].  Fix  a  small  time  increment  e  >  0  and  suppose 
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0  <  s  <t.  Utilizing  the  Markov  property  of  Z,  we  obtain 
Vij(x,s,t  +  e)  =  P(X(f  +  e)  <  x,  Z{t  +  e)  =  j\Z(s)  =  i) 

i 

=  Y  p  W  +  e)  <  x,  Z(t  +  e)  =  j\Z(t )  =  k,  Z(s)  =  i )  P(Z(t)  =  k ) 


k= 1 

e 


J]P(Z(t  +  e)  =  j\Z{t)  =  k,Z(s)  =i)F(X(t  +  e)  <  x\Z(t)  =  k,Z(s)  =  i)P(Z(i)  =  fc) 

fc=i 

t 

J^P(Z(t  +  e)  =  j|Z(t)  =  k)F(X(t  +  e)  <  x,Z(t )  =  jfe|Z(a)  =  t) 

fe=i 

i 

Y,Pkj(t,t  +  e)  P(X(t +  e)  <  x,Z(t)  =  k\Z(s)  =  i ) 

fc=i 

e 

Y/Pkjjt;  t  +  e)  Vik(x  -  r(k)e,  s,t ).  (4) 


fc=i 


Substituting  the  equality, 

{i +  ^(t)e  +  o(e),  k  =  j , 
in  equation  (4)  and  rearranging  terms,  we  obtain 


Pkj(t,t  +  £ )  = 


Vij(x,s,t  +  e)  =  [1  +  qjj(t)£  +  o(e)\Vij(x  -  r(j)e,s,t)  +  Y  Pkj(t,t  +  s)Vik(x  -  r(k)e,  s,t) 

kes\{j} 

e 

=  Vij(x  -  r(j)e,  s,  t)  +  eY<lkj(t)Vik(x  -  r(k)e,s,t )  +  o(e). 


k= 1 


Therefore, 

i 

Vij(x,s,t  +  e)  -  Vij(x  -  r(j)e,s,t )  =  e  ^  qkj(t)Vik(x  -  r(k)e,s,t )  +  o(e).  (5) 

fc=i 

Adding  and  subtracting  Vij(x,s,t )  to  the  left-hand  side  of  equation  (5)  gives 

[Uy(x,s,t  -be)  -  Uj(x,s,t)]  +  s,t)  -  Vij(x  -  r(j)e,s,t )] 

t 

=  sYQkj(t)Vjk{x  -  r(k)e,s,t )  +  o(e).  (6) 


fc=i 

Now,  dividing  each  term  of  (6)  by  e  and  allowing  e  j  0,  we  obtain 
dVij(x,s,t )  dVij(x,s,t ) 

- +  r(j)  ax 

k 


dt 


Y^kjjt)  Vik(x,  S,  t), 


(7) 


which  can  be  written  in  matrix  form  as 

dV(x,s,t)  dV(x,s,t) 


dt 


+ 


dx 


Tid  =  V(x,s,t)Q(t). 
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We  pause  here  to  remark  that  the  nonhomogeneous  equation  (3)  extends  the  Chapman-Kolmogorov- 
type  equation  derived  in  [13,  14]  to  scenarios  in  which  the  dynamics  of  the  environment’s  evolution 
depend  explicitly  on  time.  While  equation  (3)  is  more  general,  its  parametrization  can  be  problem¬ 
atic  in  practice  if  the  environment  cannot  be  continuously  observed.  Nevertheless,  the  model  can 
be  used  to  capture  these  complex  dynamics,  particularly  if  the  manufacturing  equipment  is  heavily 
influenced  by  the  effects  of  aging. 

Note  that  equation  (7),  when  expanded  term-by-term,  forms  a  system  of  partial  differential 
equations  that  must  be  solved  simultaneously.  For  example,  if  Z  has  the  state  space  S  =  {1,2} 

(i.e. ,  there  are  two  distinct  environment  states),  the  resulting  system  has  four  equations  and  four  un¬ 
known  functions  (namely  V\  \ ,  V\2,  V2\,  and  V22)  with  the  obvious  boundary  condition  Vij( 0,  s,  t)  =  0 
and  initial  conditions  Vu(x,  s,  s)  =  1  and  Vij(x,  s,  s )  =  0,  j  7^  i.  This  set  of  partial  differential  equa¬ 
tions  is  given  by 


dVn 

dt 

dV12 

dt 

dV21 

dt 

dV22 

dt 


+  r(l) 
+  r(2) 
+  r(l) 
+  r(  2) 


dVn 
dx 
dV 12 
dx 
dV 21 
dx 
dV 22 
dx 


qn{t)Vn  +  q2i(t)Vi2, 
qi2(t)Vu  +  q22(t)Vi2, 
qn{t)V2i  +  q2i(t)V22, 
q\2{t)V2l  +  <722  (i)  V22  - 


In  general,  if  the  environment  has  t  {t  <  00)  states,  then  we  must  solve  C2  partial  differential  equa¬ 
tions  with  l'2  unknown  functions.  Since  the  spatial  and  temporal  variables  (x  and  t,  respectively) 
are  weakly  coupled,  the  method  of  separation  of  variables  cannot  be  used  to  obtain  V(x,s,t )  in 
closed-form.  Therefore,  our  aim  is  to  characterize  V(x,s,t )  via  a  (Laplace)  transform  solution. 

To  that  end,  let  u  be  a  transform  variable  with  Re(u)  >  0  and  let  0  <  s  <  t.  Define  the  l  x  t 
matrix 

/•OO 

V*(u,s,t)=  /  e~uxV(x,  s,t)dx,  (8) 

Jo 

the  Laplace  transform  (LT)  of  the  matrix  V(x,s,t )  with  respect  to  the  degradation  threshold  x. 
For  later  use,  we  also  define 

V(u,s,t)=  /  e~uxV(dx,s,t), 

Jo 

the  Laplace-Stieltjes  transform  (LST)  of  V(x,s,t )  with  respect  to  x  and  note  that 


V(u,s,t )  =  uV*(u,  s,  t). 


(9) 


Now,  taking  the  LT  of  both  sides  of  equation  (3)  with  respect  to  x  and  using  the  well-known  identity 


£ 


=  uV*(u,s,t)  -  V(0,s,t), 


we  obtain 


dV*(u,  s,t) 


+  (uV*(u,  s,t)  -  V(0,s,t))Rd  =  V*(u,  s,  t)Q(t). 


dt 


(10) 


Note  that  s  is  a  fixed  initial  time  and  u  can  be  viewed  as  a  constant.  Therefore,  by  applying  the 
boundary  condition,  V(0,s,t)  =  0  (the  zero  matrix),  and  rearranging  the  terms  of  (10),  we  obtain 
a  linear  first-order  differential  equation  in  t, 


dV*(it,  s,  t) 
d  t 


V*(u,  s,  t)A(u,  t ) 


(11) 


where  A(u,t )  =  Q(t)  —  uHd,  t  >  0.  Our  objective  is  to  solve  the  following  initial  value  problem 
(IVP)  in  the  Laplace  transform  space: 


dV*(a,g,t) 

dt 

V*(u,s,s) 


V*(u,  s,  t)A(u,  t ), 
u~ll 


where  (12b)  follows  from  the  initial  condition 


Vij(x,s,s) 


1,  3  =  i 
0,  j  +  i. 


(12a) 

(12b) 


In  what  follows,  exp(A)  denotes  matrix  exponentiation  of  a  square  matrix  A,  and  the  commutator 
of  two  square  matrices  A  and  B  is  defined  to  be 


[A,B]  =  AB-  BA. 

Theorem  2  characterizes  the  solution  of  IVP  (12)  using  a  result  due  to  Magnus  [25]. 


Theorem  2  For  any  fixed  initial  time  s,  the  solution  to  the  initial  value  problem  (12)  is  the  matrix 
V*  (u,  s,  t )  given  by 

V*(u,  s,  t )  =  u~l  exp  [f l(u,  s,  i)] 


where 


and 


Cli(u,s,t) 


fl2{u,s,t) 


n3(u,s,t) 


Q(u,s,t)  =  ^2  Qk(u,s,t ) 
k= 1 


(  A(u,t i)dti 
J  S 

7}  j  [A(u,  ti),A(u,  t2)]  dt2dti 

^  L  L  I  ^A^Ujtl^^A^Ujt2^A^u,t3^  +  ^Ujt^^A^u't2^A^u,tl^dt3dt2dtl 


Proof.  Note  that  equation  (11)  can  be  easily  solved  using  an  integrating  factor  if  A(u,v) 
and  A(u,  w)  commute  for  any  pair  of  t  values  v,  w  €  M+.  In  general,  these  matrices  do  not  commute; 
however,  by  invoking  Magnus’  Theorem  (see  [3,  25]),  we  can  assert  that  solutions  of  (11)  are  matrix 
exponentials  of  the  form 

C  exp  [n(tt,  s,  t)] 
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where  C  is  a  constant  square  matrix  (with  respect  to  t),  Cl(u,  s,t )  is  the  solution  of  the  differential 
equation 


j  00  p 

Ata  =  E 

n= 0 

=  A(u,  t)  -  i  [fi,A(u,  f)]  +  [ft,  [ft, A(u,  f)]]  H - 

B„  is  the  nth  Bernoulli  number,  and  adn  is  an  operator  acting  on  a  square  matrix  defined  by  the 
commutator 

ad n(B)  =  [ft,  B]  =  ft£>  —  BCl. 

(This  notation  suppresses  the  symbols  ( u,s,t )  in  ft.)  The  series  expansion  for  the  solution  matrix 
Cl(u,s,t), 

OO 

Cl(u,s,t)  =  ^ Clk(u,s,t ), 

fc=i 

determined  by  Magnus  [25]  and  commonly  referred  to  as  the  Magnus  expansion ,  is  given  by 


Cl(u,s,t)  =  J  A(u,fi)dfi  +  i  j 


t  r  /*£i 


A(n,t2)dt2,  A(n,  ti) 


dfi 


J  J  {[A(u,t1),[A(u,t2),A(u,t3)]\  +  [A(u,t3),[A(u,t2),A(u,ti)]\)dt3dt2dti  + 


which  is  known  to  commute  with  j^Cl(u,  s,t).  Consequently,  we  may  write 


—  exp  [Cl(u,  s,  t)]  =  exp  [ft  (u,  s,  t)]  A(u,  t), 
so  the  solution  to  the  initial  value  problem  (12)  is  given  by 

V*(u,  s,  t )  =  rt_1I  exp  [Cl(u,  s,  t)]  =  u ~1  exp  [ft(it,  s,  f)]  . 


By  noting  that  V(u,s,t )  =  uV*(u,s,t),  Theorem  2  allows  us  to  write  LST  expressions  for  the 
unconditional  and  conditional  c.d.f.s  of  Tx(s),  respectively.  When  Z  is  a  NHCTMC,  it  is  well-known 
that  the  distribution  of  Z(s)  is  given  by  the  row  vector 

«(«)  =  [p(z 0s)  =  *)]*es  =  a(°)p(°’ s)  =  Q(°)  exP  Q(r)d^)  .  (13) 

To  simplify  notation  in  the  results  that  follow,  let  us  also  define 


e'(S)  =  [l(Z(s)  =  i)]i6ff 


a  row  vector  of  order  l  whose  elements  are  all  zero  except  for  the  zth  entry  which  is  unity  if  Z(s)  =  i, 
and  let  e  =  [1, 1, ....  I]7  be  a  column  vector  of  ones  of  order  i.  Note  that  ol  and  e!i  both  depend 
on  the  initial  time  s  because,  although  the  initial  cumulative  degradation  X (s)  =  0  a.s.,  we  must 
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know  the  distribution  of  the  state  of  the  environment  at  time  s  to  obtain  the  conditional  and 
unconditional  distributions  of  Tx(s)  (via  their  respective  transforms).  The  unconditional  c.d.f.  of 
Tx(s)  is  given  by 

£  £ 

F(x,  s,  t)  =  ¥(Tx(s)  <  t)  =  1-  EE  Vij(x,s,t)¥(Z(s)  =  i ) 

i=  1  j= 1 

=  1  —  a(s) V(x,  s,  t)e,  (14) 

and  its  Laplace-Stieltjes  transform  will  be  denoted  by 

F(u,s,t)=  /  e~ux  F(dx,  s,  t). 

Jo 

Likewise,  for  each  i  G  S,  define  the  conditional  c.d.f.  of  Tx(s),  given  that  Z(s)  =  i,  by 


Fi(x,s,t)  =  ¥(Tx(s)  <  t\Z(s)  =  i)  =  1  ~  y ^Vij(x,s,t) 

=  1  -  e-(s)Vr(x,s,t)e, 


(15) 


and  its  LST 


Fi 


no o 

(u,s,t)=  /  e~ux  Fi(dx,  s,t). 
Jo 


With  these  definitions,  we  can  now  state  the  following  proposition. 


Proposition  1  Suppose  the  unit  operates  in  a  time-nonhomogeneous  CTMC  environment,  { Z{t )  : 
t  >  0}  on  state  space  S  =  {1,2,...,  €}  with  time- dependent  infinitesimal  generator  matrix  Q(t). 
Then  the  LST  of  the  unconditional  lifetime  c.d.f.,  F(x,s,t),  is  given  by 

F(u,s,t)  =  l  —  a(s)exp[fl(u,s,t)]e,  Re(u)  >  0,  (16) 

and  the  LST  of  the  conditional  lifetime  c.d.f.,  Fi(x,s,t),  is 

Fi(u,  s,  t)  =  1  —  e[{s)  exp  [fi(u,  s,  f)]  e,  Re(rt)  >  0,  (17) 

where  ft(u,s,t )  is  given  by  the  Magnus  expansion  of  Theorem  2. 

Proof.  Equations  (16)  and  (17)  are  obtained  directly  by  taking,  respectively,  the  LSTs  of  (14) 
and  (15)  with  respect  to  x  and  substituting  the  result  V(u,s,t )  =  exp  [0(n,  s,t)].  m 

The  results  of  Theorem  2  and  Proposition  1  serve  to  generalize  known  results  for  homogeneous 
environments  to  a  much  broader  class  of  temporally  nonhomogeneous  environments.  Indeed,  to 
ensure  that  the  solution  matrix  V*(u,  s,  t )  is  of  the  exponential  form  exp(f2(t(,  s,  t))  we  require  only 
that  the  entries  of  Q(t)  are  Riemann  integrable  functions  of  t  that  satisfy 

(Q(r)  -  itRcf) ||2  dr  <  7 r 
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(see  Moan  [31])  where  ||  •  ||p  denotes  the  usual  p-norm  of  a  matrix.  Magnus  [25]  showed  that  the 
series  J2kLi  ^k(u,  s,  t )  will  terminate  at  a  finite  index  M  if 

=  0, 

for  all  values  of  t.  That  is,  if  A(u,  t)  and  its  integral  commute  for  all  t,  then  there  exists  a  finite 
index  M  such  that 

M 

fl(u,s,t)  =  ynk(u,a,t). 

k= 1 

For  the  special  case  when  Q(t )  =  Q  for  all  t  >  0,  we  can  state  the  following  corollary  to  Proposition 
1  which  confirms  prior  results  for  homogeneous  Markovian  environments. 

Corollary  1  Suppose  {Z(t)  :  t  >  0}  is  temporally  homogeneous,  i.e.,  Q(t)  =  Q  for  all  t  >  0.  Then 
the  LST  of  the  unconditional  lifetime  distribution  is  given  by 

F(u,  s,t)  =  1  —  a(s)  exp  [(Q  —  uRd)(t  —  s)]  e,  (18) 

and  the  conditional  version  of  the  lifetime  distribution  is  given  by 

Fi(u,s,t )  =  1  -  e'(s)exp  [(Q  -  uRd)(t  -  s)]  e.  (19) 

Proof.  If  Q(t)  =  Q  for  all  t  >  0,  then  for  u  fixed  and  any  v  £  M+,  A(u,v )  =  A(u)  = 
Q  —  uRrf.  Hence, 

[A(u,  v),  A(u,  u;)]  =  0 

for  any  pair  of  times  v,  w  £  R+.  The  terms  of  the  Magnus  expansion  of  Theorem  2  are  given  by 
fii (u,s,t)=  I  A(u,ti)dti=  f  (Q-uRd)dti  =  (Q-uRd)(t-s), 

J  S  J  S 

and  ftk{u,s,t)  =  0  for  k  >  2  since  all  commutators  in  the  subsequent  terms  result  in  the  zero 
matrix.  Consequently, 

f l(u,  s,  t )  =  s,  t)  =  (Q  —  uRd)(t  —  s). 

Equations  (18)  and  (19)  follow  directly  from  Proposition  1.  ■ 

Remark:  It  is  worth  noting  that  if  the  unit  is  placed  into  service  in  a  time-homogeneous  environ¬ 
ment  at  time  s  =  0,  then  we  can  replace  ol(s)  and  e'(s)  by  a  =  a(0)  and  e'  =  e'(0),  respectively, 
and  f l(u,s,t)  by  Q(u,t)  =  (Q  —  uRd)t.  The  corresponding  result,  expressed  as  a  Laplace  trans¬ 
form,  is  derived  by  Kulkarni,  et  al.  [21].  However,  the  methods  used  to  obtain  this  one-dimensional 
transform  of  the  lifetime  distribution  are  quite  different.  Moreover,  when  s  =  0  in  equations  (18) 
and  (19),  the  results  are  in  agreement  with  equations  (12)  and  (13)  of  [15]  when  the  effect  of  shocks 
are  ignored,  i.e,  when  A  =  0. 

When  Z  is  a  NHCTMC,  it  is  not  immediately  clear  how  to  obtain  the  nth  moment  of  Tx(s) 
(or  its  LST)  using  equations  (16)  and  (17)  because  the  transforms  and  their  inverse  transforms  are 
difficult  to  obtain  (even  through  numerical  inversion  techniques).  Suppose  we  let  F(x,s,t)  be  the 


A(u,t),  /  A(u,  r)dr 


12 


approximate  lifetime  c.d.f.  obtained  via  numerical  inversion  of  (16).  Then  by  applying  the  method 
of  integrating  the  tail,  the  approximate  nth  moment  of  Tx(s)  is  given  by 

/OO 

ntn~1[l-  F(x,s,t)]dt. 

An  alternative  approach,  demonstrated  in  [14,  15],  is  to  derive  the  LST  of  F(u,s,t )  with  respect 
to  t  and  to  evaluate  its  nth-order  derivative  at  zero.  The  resulting  expression  is  the  LST  of  the 
nth  moment  which  must  subsequently  be  inverted  numerically.  However,  a  closed-form  expression 
for  this  LST  does  not  appear  to  be  attainable  unless  A(u,  v )  and  A(u ,  w)  commute  for  any  pair  of 
times  v  and  w.  As  shown  in  this  section,  the  commutative  property  holds  when  the  environment 
is  temporally  homogeneous;  however,  it  does  not  hold  in  general.  Further  analysis  is  needed  to 
examine  the  moments  of  the  lifetime  distribution  when  Q{t )  is  not  constant. 

In  the  next  subsection,  we  illustrate  via  a  small  numerical  example  how  the  distribution  of 
Tx(s)  can  be  approximated  using  equation  (16). 


3.2  A  Numerical  Illustration 


In  general  it  is  a  difficult  task  to  compute  or  invert  the  Laplace-Stieltjes  transforms  of  Proposi¬ 
tion  1  due  to  the  complexity  of  resolving  the  individual  terms  of  the  Magnus  expansion  of  Theorem 
2.  However,  approximate  LSTs  (and  their  inverse  transforms)  can  be  obtained  by  considering  only 
the  first  term  of  the  Magnus  expansion  s,t))-  It  is  important  to  note  that  the  first  term 

cannot  give  the  whole  solution  of  (12).  Indeed,  the  remaining  terms  in  the  expansion  provide  the 
correction  needed  to  maintain  an  exponential  form  of  the  solution  as  noted  by  Blanes  et  al.  [3]. 
Nonetheless,  we  set  out  to  illustrate  the  form  of  such  an  approximation  via  a  specific  numerical 
example.  This  will  be  achieved  using  numerical  Laplace  transform  inversion. 

Suppose  the  environment  Z  is  a  NHCTMC  with  state  space  S  =  {1,2},  time-dependent  gener¬ 
ator  matrix 


Q(t) 


— 3t  3  t 
2  r  —2  r 


and  degradation  rate  matrix 


Rd  — 


1  o 

0  4 


Let  us  define  G(u,s,t )  as  the  approximate  LST  of  Tx(s)  obtained  by  evaluating  F(u,s,t)  using 
only  the  first  term  of  the  Magnus  expansion,  i.e., 


G(u,  s,t)  =  1  — 


a(s)  exp  [I2i(w,  s,  t)]  e 


where 


fii (u,s,t)=  I  A(u,  r)dr  =  f  (Q(t)  —  uR^dr. 

J  S  J  S 

It  is  not  difficult  to  show  that  for  0  <  s  <  t, 


—3  r  —  u  3  r 

2  r  —  2t  —  4  u 


dr  = 


“(f  t  +  \s  +  u)  |  (t  +  s) 

t  +  s  —  (£  +  s  +  4  u) 


(t  ~  s )■ 
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Therefore,  the  approximate  LST  of  F(x,s,t ),  with  respect  to  x,  is 


G(u,  s,  t)  =  1  —  a(s)  exp 


-( 


ft  +  fs  +  u) 
t  -\-  s 


f  (t  +  s) 

—  (t  +  s  +  Au) 


( t-8 ) 


We  arbitrarily  chose  a  fixed  initial  time  s  =  1.0.  To  compute  the  vector  o:(1.0),  we  assume  that 
Z( 0)  =  1  with  probability  0.25  and  Z(0)  =  2  with  probability  0.75  so  that  a(0)  =  [0.25,0.75]. 
Applying  (13)  we  obtain 


a(1.0) 


ck(0)  exp 


Q(r)dT 


[0.3877,0.6123], 


We  computed  the  inverse  Laplace  transform  with  respect  to  u  given  by 


G(x,s,t)  =  CUY  (u  1G(u,s,t)^j 


using  the  numerical  inversion  algorithm  due  to  Abate  and  Whitt  [1]  for  a  range  of  t  values  (t  >  s  = 
1.0)  when  x  =  2.0  units  of  degradation.  The  matrix  computations  and  numerical  inversion  algorithm 
were  coded  in  the  MATLAB  computing  environment  and  executed  on  a  personal  computer. 

For  interpretive  purposes,  we  plot  the  approximate  survivor  function  P(Ta,(1.0)  >  t)  »  1  - 
G(x,  1.0,  t)  which  represents  the  probability  that  the  unit’s  lifetime  exceeds  t,  given  that  it  was 
placed  into  service  at  time  1.0  with  no  accumulated  degradation.  Figure  1  depicts  the  behavior  of 
the  survivor  function  which  is  bounded  above  by  1,  bounded  below  by  0  and  monotone  decreasing 
in  t,  as  expected. 
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1.0 


In  the  next  section,  we  consider  an  extension  to  the  case  of  a  semi-Markov  environment  whose 
evolution  is  time- homogeneous. 

4  Homogeneous  Semi-Markov  Environment 

In  this  section,  we  consider  the  case  when  Z  =  {Z(t)  :  t  >  0}  is  a  temporally  homogeneous 
semi-Markov  process  (SMP)  with  a  Markov  renewal  process  embedded  at  transition  times  (jump 
epochs).  Let  us  define  rn  as  the  nth  jump  epoch  and  let  =  Z(t+  )  be  the  state  of  the  environment 
just  after  the  nth  transition.  The  bivariate  sequence  of  random  variables,  {(£n,rn)  :  n  >  0}  is  a 
Markov  renewal  process  with  semi-Markov  kernel  matrix,  G(t)  =  [Glj(t)\  given  by 

Gij(t)  =  P(£n_|_i  =  j,  Tn_)_ i  —  Tn  Z  t\Tn,  =  i, . . . ,  To,  Co) 

=  I®(Cn+l  =  3,  Tn+ 1  -Tn<  f|Cn  =  *)>  hJ  ^  S,  t>  0,  (20) 

where  to  =  0  and  €  S  w.p.  1  for  each  n  >  0.  Because  we  assume  the  SMP  is  time-homogeneous, 
equation  (20)  is  equivalently 

Gij(t)  =  P(Cn+l  =  j,  Tn+ 1  -  Tn  <  f  |Cn  =  *)  =  P(Cl  =  3,  U  <  *|Co  =  *) 

for  each  n  >  0.  We  assume  that  Gij(t)  is  absolutely  continuous  and,  therefore,  possesses  a  proba¬ 
bility  density  function  (p.d.f.).  The  embedded  Markov  chain,  {£n  :  n  >  0},  has  one-step  transition 
probabilities 

Pij  =  P(Cn+l  =  j|Cn  =  *)i  b  J  ^  *5*1 

and  we  let  P  =  [p^]  denote  its  transition  probability  matrix.  The  matrix  P  can  be  obtained  by 

P  =  lim  G(t), 

t— XX) 
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i.e.,  pij  =  Gij{ oo)  <  1. 

The  environment’s  holding  time  (or  sojourn  time)  in  state  i  £  S  can  also  be  obtained  directly 
from  the  semi-Markov  kernel.  Specifically,  if  the  c.d.f.  of  the  holding  time  in  state  i  is  Hi,  then 

Hi{t)  =  P(n  <  t|£0  =  i)  =  Y  G  si 

j£S 

and  its  mean  is  given  by 

/»oo 

Hi  =  xH{(dx),  i  £  S. 

Jo 

For  a  SMP,  Gij{t )  =  pijHi(t).  Assuming  P  is  irreducible  and  recurrent,  it  has  an  invariant  prob¬ 
ability  vector  7T  =  [7 Tj].  For  our  purposes,  we  assume  the  holding  time  in  each  state  is  absolutely 
continuous  so  that  the  p.d.f.  associated  with  Hi  is  denoted  by  hi,  and  we  assume  that  Hi  <  00  for 
each  i  £  S. 

As  before,  we  let  X(t)  denote  the  cumulative  degradation  up  to  time  t,  and  let  N(t)  denote  the 
total  number  of  transitions  of  Z  up  to  time  t.  Define  the  set  of  integers,  M  =  {n  >  0  :  rn+ 1  <  t}. 
Then  we  can  express  X  (t)  by 

X(t)  =  r^n){Tn+i-Tn)  +  r{Z{t))'y{t)  (21) 

n£j\T 

where  7 (t)  =  t  —  rjv(t)  is  the  age  process.  Thus,  we  can  view  {X(t)  :  t  >  0}  as  a  semi-Markov 
reward  process  with  positive,  linear  reward  function  r(-).  While  many  authors  have  analyzed 
the  asymptotic  behavior  of  X(t)  (cf.  [7,  8,  9,  17,  18]),  the  transient  analysis  of  its  distribution  is 
notoriously  difficult.  However,  by  again  exploiting  the  dual  relationship  between  X(t)  and  Tx(s),  we 
provide  a  matrix  Laplace  transform  expression  for  the  unconditional  transient  lifetime  of  the  unit. 
Because  the  SMP  environment  is  temporally  homogeneous,  we  assume  (without  loss  of  generality) 
that  the  unit  enters  service  at  time  zero  and  drop  the  dependence  of  Tx  on  s.  Likewise,  we  write 
F(x,t)  and  F)(x,t)  for  the  conditional  and  unconditional  distribution  functions  of  Tx. 

4.1  Laplace-Stieltjes  Transform  of  Tx 

Our  approach  to  deriving  the  transient  lifetime  distribution,  F(x,  t ),  is  in  the  spirit  of  techniques 
employed  by  Kulkarni  et  al.  in  [22].  Specifically,  we  will  use  a  Markov  renewal  argument  to  derive 
the  Laplace-Stieltjes  transform  (LST)  of  Tx. 

For  any  t  >  0,  whenever  Z{t)  =  j  G  S,  the  rate  of  degradation  of  the  system  is  a  positive 
number  r(Z(t))  =  r(j).  Define  a  column  vector  of  these  t  degradation  rates  by 

r  =  (r(l),r(2),...  . 

For  each  i,j  £  S  and  t  >  0  let  ipijit)  =  Gij(.t),  and  define  the  square  matrix  Let 

the  conditional  lifetime  distribution  be 

Fi(x,t )  =  P(TX  <  t|Z(0)  =  i)  =  1  —  P(A(f)  <  x\Z(0)  =  i ), 

and  let  its  LST,  with  respect  to  t,  be  denoted  by 

Fi(x,  s)  =  E  (e~sTx \Z(0)  =  1s)  =  /  e~stFi(x,dt) 

Jo 
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F(x,  s)  =  Y:  Fi(X :  s)P(Z(0)  =  i) 

2=1 

is  the  LST  of  the  unconditional  distribution  function  F(x,t)  =  P (Tx  <  t)  with  respect  to  t.  (Here, 
s  is  a  complex  transform  variable  with  Re(s)  >  0.)  Next,  we  define  the  Laplace  transform  (LT)  of 
F(x,  s )  with  respect  to  x  by 


e  uxF(x,s)dx,  Re(u)  >  0. 

To  simplify  notation  in  the  results  that  follow,  let  us  also  define  the  LT  of  ifij(t)  by 

/*oo 

Re(tu)  >  0, 

Jo 

and  for  each  i  £  S,  let 

/»oo 

H*(u)  =  /  e-^Hi(t)dt 

Jo 

denote  the  LT  of  the  sojourn  time  distribution  Hi.  Theorem  3  provides  a  closed-form  result  for  the 
(double)  Laplace  transform  F*(u,  s ).  Unfortunately,  a  one-dimensional  transform  result  appears  to 
be  available  only  if  Z  is  a  time-homogeneous  CTMC. 


Theorem  3  If  the  environment  process  Z  is  a  time-homogeneous  semi-Markov  process  with  finite 
state  space  S  =  {1,2, ...  ,£}  and  initial  distribution  a,  then 

F*(u,  s)  =  a  [I  -  */>*(«,  s)]-1[A*(u,  s)  -  H  *(u,  s)]r  (22) 

where  the  matrix  xf* (u,  s)  =  [ip*j(u,s)\  has  elements 

ij(u ,  s )  =  pij  h*(s  +  r(i)u),  i,j  €  S, 

A *(u,s)  =  diag[l/(s  +  r(l)u),  l/(s  +  r(2)u), . . . ,  l/(s  +  r(£)u)\,  and  H*(u,s)  =  diag[iL*(s + 
r(l)u),  H$(s  +  r(2)u),  ...,Hf(s  +  r(i)u)\. 


Proof.  The  result  will  be  proved  using  a  Markov  renewal  argument  by  conditioning  on 
the  first  transition  time  t\  and  the  initial  state  of  the  environment  Z( 0).  To  this  end,  note  that 


E(e“sT"|ri  =  h,Z(  0)  =  i)  = 


e  sx/r(*))  if  h  >  x/r(i ), 

e~shYlej=iPijFj(x  ~  r(i)h,s),  if  h  <  x/r(i). 


First,  we  uncondition  with  respect  to  the  first  sojourn  time  to  obtain 

~  r°° 

Fi(x,s )  =  /  E(e_sTa:|ri  =  h,Z(  0)  =  i)dHi(h) 

Jo 


[r(i)  e-shFj(x  -  r(i)h,  s)dHl(h)  +  e~sx^  [1  -  H,  (x/r(i))}  .  (23) 

i=i  Jo 
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Now,  taking  the  Laplace  transform  of  both  sides  of  equation  (23)  with  respect  to  x  and  using  the 
change  of  variable,  y  =  r(i)h  and  d y  =  r(i)dh,  we  obtain 


F*(u,s)  =  f  e  uxFj.(x,s)dx 

Jo 


£  /*oo  px  /  \ 

Tq  /  e~UX  /  s)e~sy/r{-l)hi  (  -J-  )  dydx  + 

J o  Jo  WO / 


KO 


s  +  r(z)u 

1*00 

-  /  e-uxe-S3;/,'WiLl(x/r(i))da:.  (24) 

Jo 

Using  standard  Laplace  transform  tables  (cf.  Oberhettinger  and  Badii  [33]),  it  is  known  that 


H*{r(i)u  —  s)  =  f  e 
Jo 


r(i ) 


-da;. 


(25) 


Therefore,  the  right-most  integral  of  equation  (24)  can  be  simplified  to  obtain 


F*(u,s)  = 


+ 


r(i) 


-  r{i)H*{s  +  r(i)u). 


s  +  r(i)u 

To  resolve  the  double  integral  of  (26),  we  note  that  for  some  function  g(y,  s), 


(26) 


PX  ~ 

M(x,s)=  /  Fj  (x  -  y,  s)  g(y,  s)dy 
Jo 


is  the  convolution  of  the  functions  Fj(-,  s )  and  g(y ,  s).  It  is  well  known  that  the  Laplace  transform 
of  M(x,s),  denoted  by  M*(u,s),  is  given  by 

/OO  _____ 

e~uxM(x,  s)dx  =  Fj  ( u ,  s)  g*{u ,  s), 

where  F*(u,s)  is  the  LT  of  F*(x,s)  with  respect  to  x,  and  g*(u,s)  is  the  LT  of  g(y,s).  From  (26),  we 
see  that  g(y,  s )  =  exp[— sy /r(i)\hi(y /r(i)),  and  its  LT  with  respect  to  y  is  g*(u ,  s )  =  r(i)h*(s+r(i)u). 
Hence,  equation  (26)  can  be  reduced  to 

F*{u,s)  =  ^2^-F*(u,s)r(i)h*(s  +  r(i)u)  +  ^ - r(i)H*(s  +  r(i)u) 

r(i)  J  s  +  rn)u 

j=i 

=  ^JPijFj{u,s)h*{s  +  r{i)u)+  -  r{i)H*(s  +  r{i)u).  (27) 

S  T I  i )  Ui 

3  = 1 

Now,  rearranging  the  terms  of  (27),  we  have 

1 


F*{u,s)[  1  -  pah*(s  +  r(i)u )]  -^pijF* (u,  s)h*(s  +  r(i)u)  =  r{i) 


s  +  r(i)u 


~  H*(s  +  r(i)u) 


.(28) 
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Letting  A *(u,  s )  =  diag[l/(s  +  r(l)tt),  1  /(s  +  r(2)u), . . . ,  1  /(s  +  r{£)u)]  and  H*(n,  s)  =  diag[iL^(s  + 
r(l)u),  H%(s  +  r(2)u), . . . ,  Hp(s  +  r(£)u)],  and  defining  the  ^-dimensional  column  vector,  F*(u,s)  = 
(  F*  ( u,s ) ,  P2*  (u,  s) , . . . ,  F%  (u,  s )  J  ,  equation  (28)  can  be  expressed  in  the  matrix  form 


[I  -  P  0  h *{u,  s)]F*(u,  s )  =  [A *(u,  s)  -  H*(u,  s)]r. 
where  I  is  the  identity  matrix, 


(29) 


h*(u,  s )  = 


hl(s  +  r(l)u)  hl(s  +  r(l)u )  •••  h^(s  +  r-(l)it) 
/i2(s  +  r(2)rt)  h2(s  +  r(2)n)  •••  h2{s  +  r(2)u) 


h}(s  +  r(£)u)  h}(s  +  r(£)u )  •••  h}{s  +  r(£)u) 

and  P  0  h*(u,  s )  is  the  Hadamard  product  of  P  and  h*(u,  s )  given  by 


P  ©  h*(u,  s) 


puhl(s  +  r(l)u)  p12h\(s  +  r(l)u) 
P2ih2{s  +  r(2)u)  P22h2(s  +  r(2)u) 


PuK(s  +  r{l)u) 
P2ih*2{s  +  r{2)u) 


Plih*t{s  +  r(£)u)  p£2hf(s  +  r{£)u) 


Puh}{s  +  r{£)u) 


Because  Re(u)  >  0  and  Re(s)  >  0,  and  h*(s  +  r{i)u)  =  Hj{s  +  r(i)u )  <  1  for  each  i  £  S, 
I  —  P  ©  h*(u, s)  is  invertible.  To  obtain  the  double  transform,  F*(u,s),  we  uncondition  with 
respect  to  the  initial  distribution  of  the  environment 


F*(u,  s)  =  a[I-PQ  h*(u,  s)]_1[A*(u,  s)  -  H*(u,  s)]r. 


(30) 


Now,  since  ipij(t )  =  G1- (t)  =  pVJ H'(t)  =  pijhi(t)  for  each  i,j  <E  S  and  t  >  0  in  an  SMP,  the  (i,  j)th 
element  of  P  0  h *(u,  s )  is  given  by 


+  r(i)u)  =  pij  h*(s  +  r(i)u). 

Therefore,  we  obtain  the  final  result 

F*(u,  s)  =  a  [I  -  V>*(ib  s)]-1[A*(u,  s)  -  H  *(u,  s)]r  (31) 

where  tp*(u,  s )  =  [^*-(s  +  r(i)u)].  ■ 


4.2  Illustrative  Examples 

In  this  subsection,  we  illustrate  the  form  of  the  double  Laplace  transform  for  the  case  when  hold¬ 
ing  times  are  exponentially  distributed,  hyper-exponentially  distributed,  and  Erlang  distributed. 
These  illustrations  serve  to  emphasize  the  difficulty  in  computing  the  transient  lifetime  distribution 
(or  cumulative  degradation  up  to  time  t). 
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Example  1:  Exponential  Sojourn  Times 

In  case  all  of  the  sojourn  times  are  exponentially  distributed  random  variables  with  respective 
rates  qi,  i  £  S,  Z  is  a  time- homogeneous  CTMC  on  S  with  infinitesimal  generator  matrix  Q  =  [qij\. 
The  sojourn  time  in  state  i  has  c.d.f. 


Hi(t)  =  1  -  exp(-^t), 


and  its  Laplace  transform  evaluated  at  s  +  r(i)u  is 

1  1 


H*\s  +  r(i)u)  = 
Moreover,  its  p.d.f.  is  given  by 


s  +  r{i)u  s  +  r(i)u  +  qi 
hi(t)  =  qi  exp (-qit) 


(32) 


which  shows  that 


h*  (s  +  r(i)u)  = 


Qi 


(33) 


(s  +  r{i)u  +  qi) 

It  is  worth  mentioning  that  the  transition  probability  matrix  of  the  embedded  chain,  {£„  :  n  >  0} 
is  given  by 

Qij / Qit  if  j  7^  L  Qi  7^  0) 

0,  if  j  =  i 


Pij  = 


where  qi  =  Yhj^iQij-  Substituting  P  and  equations  (32)  and  (33)  into  equation  (31),  we  obtain 


and 


I  -tp*(u,s)  = 


—q  12 

s+r(l)u+gi 

1 

_  -qe  2 

s+r(i)u+qt  s+r(i)u+qn 


-q  U 


-921 

s+r(2)u+92 


-qei 


A* (it,  s)  —  H*(u,  s)  = 


1 

s+r(l)it+9i 

0 

0 


0 

1 

s+r(2)u+q2 

0 


s+r(l)u+9i 

—q2i 

s+r(2)u+q2 


0 

0 


s+r(e)u+qe  . 

Using  (34)  and  (35),  the  elements  of  the  vector  F*(u,  s)  are  given  by 

r(i)  ,  (PJ 


F*(u,s)  = 


+  E 


s  +  r(i)u  +  qi  ,  s  +  r(i)u  +  3 

3£S\{i) 


Fj  (u,  s). 


(34) 


(35) 


(36) 


It  is  not  difficult  to  show  that  this  double-transform  result  is  equivalent  to  the  result  given  for  the 
case  of  a  continuous-time  Markov  chain  in  [21]. 
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Example  2:  Hyper-exponential  Sojourn  Times 


A  hyper-exponential  random  variable  is  a  mixture  of  exponential  random  variables  (see  for 
example  [41]).  Suppose  that  the  sojourn  time  in  state  i  6  S'  is  hyper-exponential  with  n*  phases 
such  that  the  jth  phase  has  rate  parameter  A ij,  j  =  1,2 Moreover,  let  G  (0,1)  and 
Ylj  aij  =  1  f°r  each  i  G  S.  For  t>  0,  the  c.d.f.  of  this  sojourn  time  is  given  by  (cf.  [20]) 

rii 

Hj(t)  =  1  -  exp(-Xijt), 

1= 1 

which  implies  that  its  Laplace  transform  evaluated  at  s  +  r(i)u  is  given  by 


H* 


roo 

Hi 

/  exp(— (s  +  r(i)u)t) 

1  -  ^2  exp(-A ijt) 

Jo 

l=i 

d  t 


s  +  r 

1 


roc 

— — - V'  /  exp(-(s  +  r(i)u)t)aijexp(-\ijt)dt 

r(i)u  o 


rii 

E 


s  +  r(i)u  s  +  r(i)u  +  Ajj 

The  p.d.f.  of  the  sojourn  time,  for  t  >  0,  is  given  by 

rii 

hi(t )  =  y:  qjjAij  exp(-Ajjf), 
l=i 

and  its  Laplace  transform  evaluated  at  s  +  r(i)u  is 


(37) 


roc 

h*(s  +  r{i)u)  =  /  exp(— (a  +  r(i)u)t)  /  j  ajj^ij  exp(— Xjjt)dt 

Jo  j=1 


ni 

E 


O-ij  Xij 


— (  s  +  r(i)u  +  A i 


(38) 


After  substituting  (37)  and  (38)  into  (31)  and  simplifying,  the  unconditional  double  Laplace  trans¬ 
form,  F*(u,s),  is  obtained  by 

/  \  — l  -i-i 


a 


'  n  i 

E 

\3=l 


s+r(l)u+Aij 


~P2,1 


P 11  ~P  12 


(A  g2  j\2j 

2-^  s+r(2)u+\2j 


-1 


P22 


-Pit 


~P2l 


r. 


This  result  provides  sufficient  generality  to  allow  for  each  state  sojourn  time  to  assume  a  distinct 
hyper-exponential  distribution.  Because  the  hyper-exponential  belongs  to  the  larger  class  of  phase- 
type  (PH)  distributions,  it  can  approximate  a  wide  range  of  sojourn  time  distributions.  Next,  we 
consider  the  case  of  sums  of  i.i.d.  exponential  phases  (i.e. ,  Erlang  distributed  sojourn  times). 
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Example  3:  Erlang  Sojourn  Times 

Suppose  the  sojourn  time  in  state  i  G  S  can  be  represented  by  n*  phases  (n*  >  1),  each  of  which 
has  the  same  distribution,  namely  exponential  with  mean  1/A j.  In  such  a  case,  the  c.d.f.  of  the 
sojourn  time  is  given  by 


^ 

Hi(t)  =  1  -  exp(— A t>  0. 

3=0  J' 

For  an  integer  n  >  1  and  real  number  a,  we  can  apply  the  well-known  Laplace  transform  identity 

1 


tn-lea,t 


(n  —  1)!  J  (s  —  o)n’ 
to  establish  that  the  Laplace  transform  of  Hi(t)  evaluated  at  s  +  r(i)u  is 
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The  p.d.f.  of  this  Erlang(rij,  A*)  random  variable  is 

/  \  4.\rii  —  1 

hi(t)  =  Xi  exp (—A^t) 1  1  ,  t  >  0 

tj! 

so  that  its  Laplace  transform  evaluated  at  s  +  r(i)u  is 

h*(s  +  r(i)u)  = 
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+  r{i)u  +  A 

Substituting  (39)  and  (40)  into  equation  (31),  F*(u,s)  reduces  to 
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where  the  values  a,jni,  i  €  S,  are  obtained  recursively  by 
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with  a,,  i  =  1.  Obviously,  if  n*  =  1,  for  all  i  €  S',  then  the  semi-Markov  environment  reduces  to  a 
continuous-time  Markov  chain  with  qi  =  A*,  i  G  S' . 

The  illustrative  examples  in  this  section  demonstrate  the  complexity  involved  in  characterizing 
the  reliability  (or  unit  lifetime)  of  systems  that  are  subject  to  complex  environments.  However,  it 
is  worth  mentioning  that  other  approaches  can  be  used  to  approximate  the  SMP  environment  using 
a  surrogate,  time-homogeneous  CTMC  via  phase-type  sojourn  time  distributions.  This  approach 
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has  been  outlined  and  implemented  by  Kharoufeh  et  al.  [16].  Alternatively,  to  circumvent  the 
computational  difficulties  described  herein,  it  may  be  possible  to  develop  asymptotic  results  which 
lead  to  tractable  solutions  that  are  easy  to  implement.  Khorshidian  [18]  has  recently  derived  a 
Functional  Central  Limit  Theorem  (FCLT)  using  an  approximating  martingale  process.  In  the 
next  section,  we  provide  some  concluding  remarks  and  discuss  future  work  in  this  area. 

5  Conclusions 

In  this  paper  we  have  presented  transient  reliability  indices  for  manufacturing  equipment  that 
operates  in  complex  environments  that  induce  stochastic  degradation  and  ultimately  lead  to  failure 
of  the  equipment.  We  specifically  considered  time-nonhomogeneous  Markov  environments  and  time- 
homogeneous  semi-Markov  environments.  The  mathematical  frameworks  presented  here  provide  a 
great  deal  of  modeling  flexibility,  particularly  in  scenarios  in  which  the  operating  times  in  various 
states  do  not  transition  in  a  stationary  way,  or  those  times  are  distinctly  non-exponential.  We  posit 
that  these  results  can  be  very  useful  in  scenarios  where  a  direct  mapping  between  the  environmental 
conditions  and  the  degradation  can  be  effectively  modeled.  This  is  typically  done  by  employing 
specific  physics-of-failure  models  for  the  type  of  equipment  under  consideration. 

Though  the  models  are  mathematically  sound,  they  can  be  difficult  to  implement  computation¬ 
ally.  In  particular,  for  the  two-dimensional  Laplace  transform  results,  numerical  inversion  appears 
to  be  the  best  option.  Some  good  algorithms  for  doing  two-dimensional  numerical  Laplace  trans¬ 
form  inversion  are  described  in  [2,  32].  In  the  future,  it  will  be  instructive  to  consider  asymptotic 
results,  particularly  for  the  SMP  environment  model,  to  obtain  tractable  reliability  indices.  Asymp¬ 
totic  results  may  provide  an  adequate  representation  of  the  reliability  indices,  particularly  when 
the  degradation  threshold  is  very  large.  This  may,  for  example,  correspond  to  the  case  of  highly- 
reliable  equipment  that  experiences  very  little  degradation  over  time.  The  asymptotic  results  may 
provide  a  starting  point  for  applications  requiring  real-time  updating  of  remaining  lifetime  distri¬ 
butions  (as  in  a  condition-based  maintenance  environment).  Specifically,  real-time  observations  of 
the  evolution  of  the  environment,  or  the  degradation  level,  can  be  used  to  update  the  parameters 
of  the  modulating  process  and  the  state-dependent  degradation  rates.  Tracking  the  occurrence  of 
environment  transitions  over  a  (sufficiently  long)  time  interval  will  provide  reasonable  statistical 
estimates,  so  long  as  the  number  of  environment  states  is  not  too  large. 

Finally,  a  nontrivial  task  is  the  estimation  of  the  degradation  rates  r(i)  that  describe  the  evolu¬ 
tion  of  degradation  as  a  function  of  the  environment.  Obviously,  for  this  purpose,  real  degradation 
data  is  required,  as  is  the  guidance  and  experience  of  subject  matter  experts.  Our  approach  pro¬ 
vides  a  very  macroscopic  view  of  degradation  whereas  wear  and  degradation  dynamics  can  be  very 
detailed  and  complex.  It  will  be  instructive  to  consult  with  engineers  and  scientists  to  determine 
appropriate  models  for  specific  material  types  and  certain  applications. 
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